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There is an increasing theoretical and observational evidence that the external magnetic 
field of magnetars may contain a toroidal component, likely of the same order of the poloidal 
one. Such "twisted magnetospheres" are threaded by currents flowing along the closed field 
lines which can efficiently interact with soft thermal photons via resonant cyclotron scatterings 
(RCS). Actually, RCS spectral models proved quite successful in explaining the persistent 
~ 1-10 keV emission from the magnetar candidates, the soft 7-ray repeaters (SGRs) and the 
anomalous X-ray pulsars (AXPs). Moreover, it has been proposed that, in presence of highly 
relativistic electrons, the same process can give rise to the observed hard X-ray spectral tails 
extending up to ~ 200 keV. 

Spectral calculations have been restricted up to now to the case of a globally twisted 
dipolar magneto sphere, although there are indications that the twist may be confined only to a 
portion of the magneto sphere, and/or that the large scale field is more complex than a simple 
dipole. In this paper we investigate multipolar, force-free magnetospheres of ultra-magnetized 
neutron stars. We first discuss a general method to generate multipolar solutions of the Grad- 
Schliiter-Shafranov equation, and analyze in detail dipolar, quadrupolar and octupolar fields. 
The spectra and lightcurves for these multipolar, globally twisted fields are then computed 
using a Monte Carlo code and compared with those of a purely dipolar configuration. Finally 
the phase-resolved spectra and energy-dependent lightcurves obtained with a simple model of 
a locally sheared field are confronted with the INTEGRAL observations of the AXPs IRXS 
J1708-4009 and 4U 0142+61. Results support a picture in which the field in these two sources 
is not globally twisted. 
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1 INTRODUCTION 

Among isolated neutron stars (INSs), a small number of objects, 
namely the Soft Gamma Repeaters (SGRs) and the Anomalous X- 
ray Pulsars (AXPs), share a number of common, peculiar char- 
acteristics. These are the huge spin down rates, P ^ 10~^°- 
10~^^ s/s, absence of radio emission (up to now only the AXPs 
XTE J1810-197 and IE 15 47.0-5408 have been detected in the ra- 
dio, ICamilo etal.ll2007 a''^. persistent X-ray luminosities in the 
range L ~ 10^^-10^^ erg/s and rotational periods P ^ 2-12 s, 
a very narrow range compared with that of other classes of INSs. 
Though SGRs activity is higher, both groups undergo erratic X/7- 
ray bursts with peak luminosities of ^ 10^° -10^^ erg/s and typical 
durations of ~ 0.1-1 s. Moreover, SGRs exhibit much more ener- 
getic events, the so-called Giant Flares (GFs), in which energies up 
to ^ 10^^ erg are released in a timescale of one second (see, for a 
recent review. lMereghettill2008h . Thanks to INTEGRAL ISGRI and 
RXTE HEXTE, AXPs and SG Rs are now known to be also per- 
sistent hard X-ray sources fe.g. lKuiper. Hermsen & Mendezll2004l : 



iKuiper et al.ll2006l : lMereghetti et aLl[2005l : lGotz et al.112006 '). Actu- 
ally, their energy output in the ~ 20-200 keV range may amount 
to as much as 50% of the total flux emitted above ^ 1 keV. Recent, 
deeper INTEGRAL observations have shown that the hard X-ray 
emission is highly phase-dependent and probabl y results from the 
superposition of different spectral components JPen Hartog et al.l 
2008a,b). 

At variance with radio-pulsars, the persistent emission of 
SGRs/ AXPs is ~ 10-100 times higher than their rotational energy 
losses. This, together with the lack of detected stellar companions, 
indicates that the persistent emission of these sources is unlikely to 
be powered by rotation or accretion. The (dipolar) magnetic fields 
inferred from the spin-down measurements, B ^ 10^^ — 10^^ G, 
largely in excess of the quantum critical field (Bq = 4.4 x 10^^ 
G), support the idea that SGRs and AXPs ar e ultra-magnetized NSs 
(or n iagnetars; IPuncan & Thompsonlll992l : [Thompson & DuncanI 
Il993h and their (persistent and bursting) emission is sustained by 
the super-strong magnetic field. Although other scenarios, mainly 
based on accretion from a debris disc left after the supernova event 
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(e.g. lAlpaill200ll : lEksi & Alpaj l2003l : lErtan & Alpa3l2003h . may 
still be considered, the magnetar model appears capable of explain- 
ing in a simple and economical way most of the observed properties 
of SGRs and AXPs (e.g. Woods & Thompson 2006). 

In a magnetar wit h surface field ^ 10^^ G the intern al field 
can reach - 10^'^ G dXhompson & DuncanI [19931 . 1 19951) . Since 
the poloidal and t oroidal components are expected to be in rough 
equipartition (see [Thompson et aP 1200 2. and references therein), 
the huge toroidal field stresses the crust, producing a deformation 
of the surface layers. This, in turns, induces a rotation of the ex- 
ternal field lines which are anchored to the star crust and leads 
to the appearance of an external toroidal component. The prop- 
erties of such a t wisted magnetosphere have been investigated by 
[Thompson et al.l (l2002h by nieans of a niodel analogous to tha t for 
the solar magnetic field (e.g. lLow & Loulll99Ql : IWolfsonlll995h . un- 
der the assumptions of a static, dipolar, globally twisted field and 
enforcing, as in the solar models, the force-free condition. 

A feature of (non-potential) force-free fields is the presence of 
supporting currents. As first suggested by T hompson et al.l (12002i) . 
thermal photons emitted by the star surface can scatter at the cy- 
clotron resonance on the charges flowing in the magnetosphere, and 
this can d rastically alter the primary spectrum. Recent, detailed cal- 
culations ( Lvutikov & G avriil 200d: lFernandez & Thompsonl2007l : 
iNobili. Turolla & Zanell2008a) of scattering onto mildly relativistic 
electrons confirmed this picture. Typical synthetic spectra exhibit a 
high-energy tail, superimposed to a thermal bump and closely re- 
semble the (empirical) "blackbody-i-power-law" model which has 
been routinely used to describe the m agnetars quiescent eniission 
in the ^ . 5-10 keV band (see again I Woods & Thompsonll200^ : 
iMereghettil l2008l . for a summary of observational results). This 
model has been successfully appl ied to the X-ray spectra of several 
AXPs/ SGRs bv lRea etal.l (l2008l. see also .Nobili. Turolla & Zan3 
l2008ah . providing direct support to the twisted magnetosphere sce- 
nario. 

The origin of the high-e nergy tails discovered with INTE- 
GRAL is much less understood. [Thompson & Beloborodovl (l2005l) 
analyzed different mechanisms within the magnetar model, and 
suggested that the hard X-rays may be produced either by ther- 
mal bremsstrahlung in the surface layers heated by returning cur- 
rents, or by synchrotron emission from pairs created higher up in 
the magnetosphere. Quite in terestingly, tearing & Hardind fcoovb 
and lBaring & Hardi"na (l2008l) have recently proposed a further pos- 
sibility, according to which the soft gamma-rays may also origi- 
nate from resonant up-scattering of seed photons, if a population 
of h ighly relativistic electrons is pr esent in the magnetosphere (see 
also lNobih. Turolla & Zan9l2008bh . 

Since the conduction current in a sheared magnetosphere is 
oc V X the particle density and hence the optic al depth to reso- 
nant scattering depends on the shear ( Thomp son et al.l2002l) . More- 
over, the global topology of the magnetic field influences scattering 
even for fixed shear. A globally twisted dipolar field will in gen- 
eral give rise to a different spectrum from that of a globally twisted 
quadrupolar field with the same twist angle. Because the spatial dis- 
tribution of charges is not homogeneous jThompson et al.ll2002h . 
changes in shear and/or field topology are going to produce differ- 
ences not only in the spectral shape, but also in the pulse shape of 
the received radiation. Some complicated pulse profiles from SGRs 
have already been explained by invoking higher o rder multipoles 
(e.g. [Thompson & Duncan|[200ll : lFeroci et al.ll200ll) . Moreover, re- 
cent data both for AXPs and SGRs seem to poi nt towards the pres- 
ence of a localized, rath er than global, twist fe.g. lWoods et al.l2007l : 
iPerna & Gottheljl2008l) . However, no investigation which includes 



multipolar components or localized twists have been presented up 
to now. 

In this paper we discuss how globally twisted magnetostatic 
equilibria can be derived in the case of higher order, axially sym- 
metric multipolar fields. In particular, explicit (numerical) solu- 
tions for quadrupolar and octup olar fields are presented. We use the 
Monte Carlo code developed bv lNobili. Turolla & Zan3 (l2008ah to 
investigate the properties of the emerging spectrum and pulse pro- 
files for higher order multipolar fields. The paper is organized as 
follows: in 32 we introduce the globally twisted model and de- 
rive the solutions for each axially symmetric multipole; an ana- 
lytical solution, valid in the case of dipolar and quadrupolar fields 
for small shear, is presented in Appendix [A] Monte Carlo spectra 
and lightcurves obtained with different force-free magnetospheric 
configurations are discussed in 33 where a comparison with the 
timing properties of the hard X-ray emission from the AXPs IRXS 
J1708-4009 and 4U 0142-1-6 is also presented. Discussion follows 



2 GLOBALLY-TWISTED AXISYMMETRIC MODELS 

In this section we present the basic equations that we use to derive 
globally-twisted ma gnetosphe r ic mo dels. We follow clos e ly the 
approach outlined in| WolfsonI Il995h and lThompson et all (l2002h 
(hereafter W95 and TLK, respectively), who considered the global 
twist of a dipolar field. As we show below, the same approach can 
be used to compute twisted fields for multipoles of arbitrary or- 
der. As in W95 and TLK, we restrict to magnetostatic, force-free 
equilibria. In the case of a low density, static plasma, in fact, in the 
standard MHD equation 

dv 

P-^ + p{v • V)^; = - Vp + pg+j X B, 

where p and v are the plasma density and velocity, the velocity and 
gravity terms can be neglected. With the further hypothesis that the 
(plasma) pressure force is small with respect to the Lorentz force 
j X B (where j is the current density), the equation reduces to 
j xB = 0. Since we are interested in stationary configurations, the 
Ampere-Maxwell equation simplifies to V x S = (47r/c)j. From 
the two previous conditions the usual expression for the force-free 
condition is recovered 



{V X B) X B = 0. 



(1) 



Our aim is to construct an axisymmetric, force-free field by 
adding a defined amount of shear to a potential field. In accordance 
with both TLK and W95, we choose to use the flux function ^ 
to express the poloidal component of the field. Axisymmetry is en- 
forced choosing a function independent of the azimuth (we use a 
spherical coordinate system with the polar axis along the magnetic 
moment vector). Thus, the general expression for an axisymmetric 
field is: 



B = — + B^{r,0)ecty. 



(2) 



where is the unit vector in the (j) direction. 

By inserting the previous expression into the force-free condi- 
tion (equation [Tl) one obtains two independent scalar equations for 
the flux function, and the toroidal component of t he field, 3^. 
The former requires to be a function of ^ only (ILow & Loul 
ll99Qh . thus we can write 



(3) 
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where F is an arbitrary function. Introducing the previous expres- 
sion into the second scalar equation leads to the Grad-Schluter- 
Shafranov (GSS) equation 



■■ 



(4) 



(here and in the following /i = cos 0). 

The GSS equation can be reduced to an ordinary differential 
equation by making suitable assumptions on the dependence of ^ 
on t he coordinates. Fol lowing a classical approach to this problem 
(e.g. lLow&Loulll99QL W95, TLK), we assume separation of vari- 
ables and choose the flux function ^ in the form 



Rns 



(5) 



where /(/x) is a function of the colatitude 0, Rns is the stellar 
radius and = Bpoiei?Ars/2, as in TLK. The requirement that 
all the components of B have the same radial dependence implies 
that 







C 



(6) 



where, for later convenience, we have expressed the multiplicative 
constant in terms of a parameter C and of the radial exponent p. 
Recalling equation (|2]), one can explicitly write the magnetic field 
as a function of / 



pole 



Rns 



Cp f 
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' sin ^ ' V P + 1 sin 



i+i/pi 



(7) 



where a prime denotes derivation with respect to fi. Finally, taking 
into account equations ^ and the GSS equation becomes for 
the case at hand 



(i-M')r+p(p+i)/+c/ 



l + 2/p , 







(8) 



which is a second order ordinary differential equation for the an- 
gular part of the flux function. Its solution, once a suitable set 
of boundary conditions has been supplied (see § I2.1K completely 
specifies the external magnetic field. 

Besides controlling the radial decay, the parameter p also fixes 
the amount of shear of the field. In fact, recalling the definition of 
shear angle (W95, TLK) 



A(I)NS 



B(h 



C 



p{l+p) 



1/2 



1 



-d/i, 



(9) 



it is immediate to see that different values of p correspond to fields 
with different shear. Actually, as it will be discussed later on, the 
effect of decreasing p is to increase with respect to the other 
components, and consequently to increase the shear. 

As it is apparent from equation 0, the limiting case p ^ 
results in a purely radial field whose field lines are directed either 
outwards or inwards (whence the name of split monopole). The di- 
rections of the field lines divide the sphere (i.e. the star) into several 
zones, and we can imagine field lines of opposite directions con- 
necting at radial infinity (W95). Split monopole fields obtained for 
different multipolar orders differ from one another in the number 
of zones into which the sphere is split. 



Figure 1. Zonal harmonics are used to visualize multipolar field topologies. 
From top left to bottom right: multipoles with order from one to four. The 
positions of the degenerate poles are evident. 



2.1 Boundary conditions 

In order to solve equation ^ we need to provide a suitable set of 
boundary conditions. For a dipolar field there are just two poles 
and the star is divided into two hemispheres, i.e. the portions of 
the surface ^ ^ 7r/2, 7r/2 ^ ^ tt, respectively, and ^ 
^ 27r. Higher order multipoles have in addition degenerate poles, 
which are loci of constant colatitude, and correspond to circles on 
the surface (see Fig.[T]). For the quadrupole, for example, there are 
two poles on the magnetic axis and the equator is a degenerate pole. 
In the case of a generic multipolar field, the magnetic hemispheres 
of the dipole are replaced by a number of regions, each limited by 
two consecutive values of the colatitude, Oi, Oi+i, for which the 
field is purely radial (the pole) and has vanishing radial component 
(the magnetic "equator"), i.e. Bo{Oi) = and Br{Oi^i) = 0. We 
refer to these zones as to regions of unipolar Br. 

Because of the north- south symmetry of unsheared multipolar 
fields, which is assumed to hold also for their twisted counterparts, 
the integration domain is restricted to ^ /i ^ 1. In analogy with 
the "composite magnetic fields" of ILow & Loul (ll99 Q), integration 
is performed piecewise, going from one pole to the next one. In 
each interval, the boundary conditions for equation ^ are deter- 
mined by the requirements that: i) B is purely radial at each pole, 
and ii) the intensity is not modified by the shear. The latter con- 
dition can be enforced either by assigning the field strength at one 
pole, or the magnetic flux out of each region of unipolar Br ; in both 
cases the value must not change with the shear. For multipoles of 
odd order (dipole, octupole, . . . ) the first sub-domain starts at the 
geometrical equator, /x = 0, which is not a pole. In these cases the 
boundary condition at /x = reflects the N-S symmetry of the field 
and translates into f{0) = (TLK). 

The magnetic field in the entire interval (0 ^ /i ^ 1) is then 
obtained by assembling the solutions computed in the various sub- 
domains with the same value of p, to ensure that the radial depen- 
dence of B is the same at all co-latitudes. Since the GSS equation is 
a second order ODE, imposing three boundary conditions implies 
that the parameter C is an eigenvalue of the problem and depends 



4 Pavan, L., Turolla, R., Zane, S. andNobili, L. 



on the radial index p, C = C{p). Each multipolar potential term 
satisfies the GSS equation with C = and is associated to an in- 
teger radial index po; it is po = 1 for the dipole, po = 2 for the 
quadrupole and so on. In this particular case, the equation is linear 
and admits analytical solutions of the form (W95) 



/po(m)oc Vl-M'-Ppo(M), (10) 

where (/j) is the associate d Legendre function of the first kind 
(lAbramowitz & Stegunll 19721) . 

Having established a set of boundary conditions, the GSS 
equation can be solved for different values of p, building a se- 
quence of models characterized by a varying shear. As shown by 
ILow & Loul (ll99Cll) . such a sequence share the same topology, i.e. 
the same number and position of poles. The only permitted values 
of the radial index are p ^ po- 



2.2 Dipolar fields 

The generating function of a pure dipole is 



Po=i 



1 



(11) 



thus from equation ([7]) it follows that the condition at the pole /i = 
1 is /(I) = 0, to which the symmetry condition f'{0) — must be 
added. The third condition is set either specifying the field strength 
at one pole, f\l) = 1 (as in TLK), or requiring a constant flux, 
/(O) = 1 (W95). The sequence of sheared dipoles has radial index 
O^p^l. 

Sheared dipole fields have been discussed by W95 and TLK. 
These investigations, however, used a different boundary condition 
(see above) and we verified that the numerical solution of equa- 
tion ^ produces quite different results for /(/i; p) and even more 
diverse eigenvalues C (p) in the two cases (see Fig.|2]). Still, one ex- 
pects the magnetic field to be the same in both cases, since the two 
boundary conditions are physically equivalent. Actually, a direct 
comparison of the numerical solutions shows that /tlk / fw95 — 
Itlk / fw95 — constant for each p. Denoting by X(p) this con- 
stant ratio, one can write the expression of the field in the two cases 
using equation 



Bw95 OC 



f ■) • /I 

sm u 



Cw95 



P+1 



1/2 jl + l/p 



sm f 



Ct 



p+1 



(12) 

1/2 j^y+i/p' 
sinO 



where the radial dependence has been omitted. The two fields have 
the same topology if and only if 



(13) 



We checked that the previous condition is indeed satisfied by our 
numerical solutions with a relative accuracy ^ 2% for nearly all 
values of pB The eingevalues and the ratio A for different values 
of p are reported in table [T] 

The magnetic field, then, has the same topology in both cases, 
and the two solutions differ only for a multiplicative factor (which 
changes withp); i.e., it is Btlk = \Bw9b, with 1 ^ A ^ 2. The 
great diversity in the eigenvalues, especially for p ^ is in fact 

^ The error becomes larger for p ~ po because A V^wm 
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Ctlk 




Cw9b 
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0.97 


0.13 




0.11 


1.02 


0.89 


0.41 




0.45 


1.06 


0.82 


0.63 




0.79 


1.10 


0.74 


0.78 




1.14 


1.15 


0.67 


0.86 




1.51 


1.21 


0.59 


0.86 




1.91 


1.26 


0.52 


0.79 




2.35 


1.33 


0.45 


0.64 




2.87 


1.39 


0.38 


0.43 




3.52 


1.47 


0.30 


0.22 




4.41 


1.56 


0.22 


0.06 




5.80 


1.65 


0.15 


3.9 X 10- 


-3 


8.52 


1.75 


0.07 


4.3 X 10- 


-7 


17.48 


1.88 


0.02 


1.55 X 10- 


-25 


89.78 


1.97 



Table 1. The eigenvalues C(p) for the two different sets of boundary con- 
ditions and the ratio A = Jtlk / fw9b for < p < 1. 
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Figure 2. Upper panels: the function /(/i) and the eigenvalue C for differ- 
ent values of p with the prescription of constant Bpoie (the dashed line in 
the left panel is the analytical solution for p 0). Lower panels: same, but 
for constant flux. 



balanced by the different behavior of the function /(/i). While for 
constant 5poie, it is 1 ^ max / ^ 2, in the case of constant flux 
max f — 1 for any p. In the limit p ^ the proportionality of the 
two solutions for / can be recovered analytically. In fact, it can be 
shown that in the split-monopole limit the generating functions are 

/(/i) = 1 — for constant flux 

/(/i) = 2 (1 — |/i|) for constant Bpoie 

2, in agreement with the numerical result (see 



which gives A(0) 
again table [T])- 



2.3 Higher order multipoles 



Since equation ^ is the force-free condition for a generic ax- 
isymmetric field (as given by equation IS), axisymmetric globally- 
twisted multipoles can be found solving again equation (|8]), sub- 
ject to the boundary conditions discussed in ^2.11 For an untwisted 
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0.0 0.2 0.4 0.6 




Figure 3. Generating functions and eigenvalues for the quadrupolar fields. 
Different curves correspond to different values of p from p = 2 (lower 
curve) to p = 0.2 (upper curve). 




Figure 4. The shear angle A(j)Ns as a function of p. The three curves cor- 
respond to dipolar quadrupolar and octupolar fields (from left to right). 



quadrupolar field the generating function is 
/po=2 = M (1 -M^)- 



(14) 



Within our integration domain ^ /x ^ 1, the poles are located at 
11 — (degenerate) and fi — 1. The boundary conditions are then 



/(I) 
/(O) 

/(I) = 



/ 



(V3 



2 

37!' 



(15) 
(16) 

(17) 



the latter two conditions enforce either constant field strength at 
the pole or constant flux, and, as discussed in the previous sec- 
tion are equivalent. However, on a numerical ground, we found 
that for higher order multipoles the constant flux boundary con- 
dition is highly preferable and has been used in the calculations 
presented here. The field in — 1 ^ /i ^ is obtained by symme- 
try. The quadrupolar angular functions /(/x) and eigenvalues C(p) 
are shown in Fig. [3] for different values of p. The shear angle as a 
function of p, equation (|9]), is shown in Fig. |4] (middle curve). 

In the case of an untwisted octupole, instead, the generating 
function is 



/po=3 = ^(l-/i^)(5/i^ - 1) 



3 0.0 




0.0 0.2 0.4 0.6 0.8 1.0 




Figure 5. Same as Fig. [3] but for octupolar fields; here p is in the range 
[0.2,3]. 



and the poles are located at /x = 1,1/ (again restricting to the 
range ^ /x ^ 1). In order to compute the sheared field, equation 
([8]) needs to be solved in two separate intervals, ^ /x ^ 1 / 
and 1/^5 ^ /i ^ 1. Taking into account that /x = is not a pole 
in the present case, the boundary conditions are 



/'(O) 
/'(1/75) = 



= 
= 
2 

7! 



/(O) 



in the range ^ /x ^ 1 / a/5 and 

/(i) = o 

/(1/75) = 
/'(I) = -2 or / 



(19) 
(20) 

(21) 



(22) 
(23) 

(24) 



in the range 1 / a/5 ^ /x ^ 1. The numerical solutions are shown in 
Fig- Eland the dependence of the shear angle on p can be read from 
Fig. |4] (rightmost curve). 

Despite the numerical solution of equation ^ poses no prob- 
lems, having analytical expressions for the twisted field compo- 
nents may prove handy in some applications. In Appendix lAl we 
derive approximated analytical expressions for the generating func- 
tion /(/x) for dipolar and quadrupolar configurations valid for small 
shear (po — P ^ 1), compare them with the exact numerical solu- 
tions and discuss their ranges of validity. 



3 SPECTRA AND LIGHTCURVES 

In this section we use th e numerical code described in 
iNobili. Turolla & Zand (l2008al . NTZ in the following) to explore 
the effects of different sheared multipolar fields on the emergent 
spectra and lightcurves of magnetars. Although the external field of 
a magnetar is unlikely to comprise a single higher order multipole, 
investigating spectral formation when the field is one (twisted) mul- 
tipole offers the opportunity to explore the effects of magneto- 
spheric currents localized on spatial scales smaller than that implied 
by the (twisted) dipole. An example is that of a sheared, localized 
component which appears as a consequence of some form of ac- 
tivity and adds up to a global, (quasi)potential dipolar field. The 
currents responsible for the resonant scatterings are provided only 
by the sheared field. The case in which the localized component is 
modeled in terms of the polar "lobe(s)" of an octupole, is discussed 



in ^3.31 We base our model on the scenario envisaged by TLK 
(18) 

(see also iLvutikov & Gavriilll2006l : [Fernandez & Thompso nl l2007l : 
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Figure 6. Optical depth to resonant cyclotron scattering at the resonant ra- 
dius as a function of the colatitude for quadrupolar configurations; each 
curve is labelled by the value of the shear angle A(f)Ns- 



iNobili. Turolla & Zanell2008aL for more detailed calculations), ac- 
cording to which thermal photons originating at the star surface un- 
dergo repeated scatterings with the charge carriers (electrons, ions 
and possibly pairs; see also Beloborodov & Thomp son 2007) flow- 
ing along the field lines. These investigations were based on non- 
relativistic computations, and therefore necessarily restricted to the 
low-energy (< 10 keV) emission. Furthermore, they were based on 
the dipolar, globally-twisted magnetosphere of TLK0. 



3.1 Globally twisted multipoles 

In order to gain some insight on the properties in the emitted spec- 
tra when the magnetosphere is threaded by twisted, higher order 
multipoles, we plot in Figs. [6] and |7] the optical depth Tres to reso- 
nant scattering corresponding to quadrupolar and octupolar fields. 
This is given by 



-(1 + COS Oks) 
4 



1/2 jl/p 



2+p 



(25) 



where v is the charge velocity and OkB is the angle between the 
primary photon (assumed to move radially) and the magnetic field 
at the scattering radius (TLK). Since v ^ c, (v/c)Tres ^ Tres and 
we make no further distinction between these two quantities. By 
comparing these curves to those in Fig. [8] we can clearly see that 
the optical depth for higher order multipoles is sensibly different 
from that of a dipolar configuration, implying that the overall spec- 
tral properties, and in particular those of the pulse phase emission, 
may be significantly affected. We remark that the typical radius at 
which resonant scattering occurs is different for the different multi- 
poles. Assuming that the surface strength is comparable, scattering 
in higher order multipolar fields takes place closer to the star with 
respect to the dipole because the field decays faster with radius. 
To further investigate this, we calculated a number of synthetic 




Figure 7. Same as in Fig.|6]but for an octupolar field. 




^ The simplified, analytical model of iLvutikov & Gavriill l2006h did not 
assume a precise topology for the magnetic field. 



Figure 8. Same as in Fig.|6]but for a dipolar field. 



spectra and lightcurves by using the non-relativistic Monte Carlo 
code by NTZ. A comparison among spectra produced by a glob- 
ally twisted dipolar, quadrupolar and octupolar magnetosphere is 
shown in Fig.[9l Spectra have been computed for the same values of 
model parameters (blackbody temperature kT^y = 0.5 keV, temper- 
ature and bulk velocity of the magnetospheric electrons Tei — 30 
keV, I3huik — 0.5, polar magnetic field strength Bpoie — 10^^ 
G), and for the same shear angle (A^a^s = 1.2 rad, which corre- 
sponds to p = 0.8, 1.6, 2.5 for the dipole, the quadrupole and the 
octupole, respectively jj. Spectra have been computed by collect- 
ing photons over the entire observer's sky (i.e. by angle- averaging 



^ Although these values of the model parameters can be regarded as typical 
(see NTZ), the present choice has only illustrative purposes. Other combi- 
nations of the parameters are equally possible and, provided that electrons 
remain mildly relativistic, will give similar results. 
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Figure 9. Monte Carlo spectra for globally twisted dipolar (solid line), 
quadrupolar (dashed line), octupolar (dash-dotted line) force free magne- 
tospheres; the shear angle is the same in all three cases (A(/)ivs=1.2 rad). 
The solid red line is the seed blackbody spectrum. The total number of pho- 
tons is ~ 4 X 10^ in all the simulations. 

over all viewing directions). The most prominent feature in Fig. [9] 
is the higher comptonization degree induced by the quadrupolar 
(and dipolar) field with respect to the octupolar one. This can be 
understood in terms of the different spatial distribution of the scat- 
tering particles (see Figs.[6j[7l[8]) and of the different efficiency of 
scatterings, the latter depending on the (average) angle between the 
photon direction and that of the flowing currents. Upscattering is 
more efficient in regions where the currents move towards the star, 
i.e. close to the magnetic south pole(s), because collisions tend to 
occur more head-on (see also NTZ). For a dipolar field the more 
favourable situation (i.e. large optical depth and currents flowing 
towards the star) arises for /i ^ —0.3 (compare with the spectra at 
different viewing angles in Fig. 1 of NTZ; the one at Gs = 116° is 
the more comptonized). For the quadrupole returning currents are 
localized around /a ^ (the geographical equator which is a de- 
generate south pole) and there are two regions with large optical 
depth in ^ /i ^ 1 (see fig. [6]). The one at /a ^ 0.3 is closer to 
the south pole (/i = 0) than in the case of the dipole, for which 
it occurs at /i ~ —0.3 while the south pole is at /x = —1. The 
reason for which octupolar twisted fields produce less efficient up- 
scattering is that, despite there are two maxima of the optical depth 
located quite close to the south pole (at /x ~ 0.2, 0.6, the south 
pole is at /i = 1/a/5 ^ 0.45), the relatively large curvature of the 
field lines makes the angular extent of the region where currents 
are inflowing narrow. As a consequence most photons scatter with 
electrons moving at large angles and this results in steeper spectra. 

3.2 A simple localized twist model 

There is now observational evidence that, at least in some cases, the 
magnetospheric twist in both SGRs and AXPs could b e localizedjn 
restricted regions of the magnetosphere. For instance, I Woods et^ 
found a certain degree of hysteresis in the long-term evo- 
lution of SGR 1806-20 prior the emission of the giant flare in De- 
cember 2004, with a non trivial correlation between spectral and 
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Figure 10. Spectra obtained with the Monte Carlo code for a locally twisted 
(A(/)jvs = 1-5 rad, p = 2.3) octupolar force-free magnetosphere. The 
shear is equally distributed on both polar lobes, while in the equatorial 
zone the magnetosphere is potential. Different curves correspond to differ- 
ent viewing angles (blue =0°, red Ss = 60°, green Ss = 104.5°, 
black Ss = 180°; the dashed line is the seed blackbody spectrum. 

timing properties that may be interpreted if only a small bundle of 
magnetic field lines is affected by the shear. A further case is pro- 
vide d by the spectral evolution of the transient AXP XTE J1810- 
197 dp erna & Gotthelfll2008l : iBernardini et al.ll2008h which seems 
to require a twist concentrated towards the magnetic axis, giving 
rise to a polar hot region, possibly with a meridional temperature 
gradient. 

Motivated by this, we investigated the possibility to model 
a localized twist, by constructing a solution in which the shear 
changes with the magnetic colatitude of the field line foot point. 
Although the configurations presented in the previous sections are 
globally twisted, the octupolar solution, which requires piecewise 
integration in two domains (see ^2.3K can be used to construct a 
field with a non- vanishing twist at low magnetic co-latitudes. This 
is done by superimposing a sheared octupole for ^ /x ^ 1/ V5 
to a potential one in 1/ a/5 ^ /x ^ 1 and is equivalent to solve 
equation ^ in the entire domain with = p*H{fi — 1/a/5), 
where H is the Heaviside step function and p* a given constant. 
We note that, although not physically self-consistent, such a field is 
force-free (see the discussion in 

The Monte Carlo spectra for different viewing angles are 
shown in Fig.[To]in the case in which the field has equatorial sym- 
metry, i.e. the shear A(/)Ars = 1.5 rad (p = 2.3) is applied on both 
polar lobes while the equatorial zone is permeated by a potential 
octupole. The spectra for the more interesting case in which the 
same shear is confined only around one pole are plotted in Fig.fTTI 

3.3 Timing and spectral properties of magnetars high- energy 
emission 

Spectral models which account for different magnetospheric con- 
figurations hold the potential to reproduce not only the gross fea- 
tures of the observed spectra but also the subtler properties which 
are revealed by the combination of very high-quality spectral and 
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Figure 11. Same as in Fig.fTolbut for the shear localized only at upper north 
pole. 



timing data. While a complete application is beyond the scope of 
this paper, here we consider our results in the context of the spec- 
tacular phase-dependence which has been recently discovered in 
the the ha rd X-ray tails of the two A XPs IRXS J1708-4009 and 4U 
0142+61 jPen Hartog etani2008al lbl). These deep INTEGRAL ob- 
servations have shown that, in both these sources, there are several 
different pulse components (at least three) with genuinely differ- 
ent spectra. The hard X-ray spectrum gradually changes with phase 
from a soft to an hard power law, the latter being significantly de- 
tected over a phase interval covering ^ 1/3, or more, of the period. 

In order to see how these features can, at least qualitatively, 
be explained by our models, let us consider the locally twisted con- 
figurations discussed in ^ 13. 2 1 which, among those presented so far, 
provide the most significant variations of the magnetic topology 
with the colatitude. Let us introduce two angles, x 6 which 
give, respectively, the inclination of the LOS and of the magnetic 
axis with respect to the star spin axis. This allows us to take into 
account for the star rotation and hence derive pulse profiles and 
phase-resolved spectra. Because of the lack of north- south symme- 
try, it is ^ X ^ while ^ spans the interval [0, 7r/2] (see NTZ 
for further details). For each viewing geometry, and for different 
values of the shear, we can now compute the optical depth (equa- 
tion l|25l ) as a function the rotational phase 7 (0 ^ 7 ^ 27r). A 
few examples are shown in Fig.[T2]for the case in which the twist 
is localized on two polar caps. 

In connection with the spectral evolution with phase observed 
in IRXS J1708-4009 and 4U 0142-^61, the more favourable cases 
are those in which it is {v/c)rres > 1 for roughly one third of 
the period. This is because resonant scattering over a population of 
(relativistic) electrons is then expected to be most efficient in pro- 
ducing a hard tail over the right phase interval, while the decrease 
of the depth at other phases results in a softening of the spectrum. A 
complete exploration of the parameter space aimed at searching for 
all configurations for which the previous condition is met is beyond 
the purposes of this paper. Just for illustrative purpose, let us con- 
sider one of those, i.e. an octupolar field with shear A^at-s = 1-5. 
By assuming this value and taking ^ = 32° x — 140°, we then 
computed phase resolved spectra and energy dependent lightcurves 
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Figure 12. The depth (v/c)Tres versus rotational phase for an octupolar 
field with twist located on the two polar caps. Each panel refer to different 
values of the geometrical angles x and ^ (see text for details). The curves in 
each panel are for Acfy^^s = 1.22, 1.50, 1.78, 1.96 from bottom to top. 



by using our Monte Carlo code. Since we are dealing with photon 
energies ^100 keV, at which electron recoil and relativistic effects 
may become important, we performed the runs using a relativistic 
version of the code ( Nobili . Turolla & Zanel l2008b. Nobili, Turolla 
& Zane, in preparation). Results are reported in Figs. [T3l and [T5l 
Again in the spirit of probing the potentialities of our model, rather 
than a presenting a detailed fit, the parameter values in the Monte 
Carlo runs are the same as those adopted in ^3.21 B^oie — 10^^ G, 
/^fe^^fc = 0.5, Tel = 30 keV and = 0.5 keV. We remark again 
that these values are not preferential, and that the results presented 
below do not critically depend on the model parameters. 

As it can be seen from Fig. (TS] for a magnetic configuration 
with the shear concentrated in a single lobe, resonant comptoniza- 
tion gives rise to a hard tail which is quite pronounced at the peak 
of the pulse while it is depressed by almost an order of magnitude 
at pulse phases close to the minimum of the hard X-ray lightcurve. 
This is similar to w hat observed in AXP IRXS J1708-4009 and 
AXP4U0142-h61 ( Den Hartog et al.l2008allbl) . This model also pre- 
dicts a considerable variation of the pulsed fraction with energy, 
ranging from a few percent below 2 keV, to a few tens of percent 
from 2 to 10 keV and up to 90% in the harder part of the spectrum. 
The comparison between modelled and observed phase-resolved 
spectra is beyond the scope of the present investigation. 

For comparison we show in Fig.Othe same results computed 
by assuming a globally twisted dipolar magnetosphere, with the 
same shear and taking the same parameters for the viewing geome- 
try and the Monte Carlo run. It is clear that the pulse resolved spec- 
tra and lightcurves obtained in this case are completely different 
and can not reproduce those observed from the two AXPs. In this 
case, the hard part of the spectrum show very little variation with 
the rotational phase and the lightcurve dependence on the energy 
is opposite to that shown in Fig. (TS] with a larger pulsed fraction 
expected in the soft band. Spectra and lightcurves obtained with 
a local twist applied at both polar regions produce results which, 
again, are not in agreement with observations (see Fig.fTSt. 



Topology ofmagnetars external field. L 9 




Figure 13. Synthetic lightcurves (upper panel) and phase-resolved spectra 
(lower panel) for a twisted magnetosphere with shear applied only to one 
polar lobe of an octupolar field (see text for details). In the upper panel the 
solid, dash-dotted and dashed lines refer to the pulse profiles in the 0.5-2, 
2-10 and 10-100 keV bands respectively. In the lower panel curves with 
different colors give the spectra at various phases; the color code can be 
read at the bottom of the upper panel. 



4 DISCUSSION AND CONCLUSIONS 



Figure 14. Same as Fig.[T3]for a globally twisted dipolar field. 



As first discussed by iThompson et alj (l2002h . the external mag- 
netic field of a magnetar likely possesses comparable poloidal 
and toroidal components. Twisted magnetospheres around ultra- 
magnetized neutron stars have been shown to play a crucial 
role in shaping the emergent spectrum of SGRs/AXPs qui- 
escent emission through efficient resonant scattering of ther- 
mal p hotons onto the charge carriers flowing along th e field 
lines ([Lvutikov & Gavrii l 2006; Fernandez & Thompson l2007l: 
iNobili. Turolla & ZandlfoOSah . 

In this paper we tackled the problem of constructing sheared 
magnetic equilibria more general than a dipole. We have show 
how sheared multipolar fields of arbitrary order can be com- 



puted by generalizing previous results by 'Wolfson ("1995') and 
Thompson et al. ( 2002). In order to assess the effects of differ- 
ent external field topologies on the emitted spectrum and pulse 
profiles we run a number of Monte C arlo simulations, using the 
code of iNobili, Turolla & Zand (l2008al) , and compared the results 
to those of a sheared dipolar field. Not surprisingly, the overall 
spectral shape does not change in going from a dipole to higher 
order multipoles and can be always described in terms of a "black- 
body plus power-law". There are, however, quite substantial dif- 
ferences among the multipoles in the spectra viewed at different 
angles. These are mainly due to the different particle distribution 
in the magnetosphere which is directly related to the assumed field 
topology. 

The case of an octupolar field has a special interest because 
it can be used to mimic a twist localized in a region close to the 
magnetic pole(s), and hence to investigate the properties of spectra 
produced in locally twisted magnetospheres. We have computed 
model spectra and lightcurves for the cases in which the twist 
is confined to one or both polar regions (each region has semi- 
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Figure 15. Same as Fig.[T3]for the case in which the shear is applied to both 
polar regions. 



aperture ^ 60°), by assuming that only the polar lobes have 
a non- vanishing shear while the equatorial belt is potential. Quite 
interestingly, a twist confined to a single lobe is the only config- 
uration, among those we have explored, that is able to reproduce 
the main features of the high-energy 10-200 keV) emission ob- 
served with INTEGRAL from the AXPs IRXS J1708-4009 and 4U 
0142-1-61, in particular the large variation in the pulsed fraction at 
different energy bands (Den Hartog et al. "2008 alibi). 

All magnetic equilibria we discussed in this paper are glob- 
ally twisted, axially symmetric multipolar fields. Of course, these 
configurations are far from being general and, even restricting to 
axial symmetry, represent only a subset of the solutions of the 
force-free equation. The magnetic field of a magnetar is likely to 
be quite complex. Modelling it in terms of single multipolar com- 
ponents offers a way of gaining insight on the general properties of 
the magnetosphere but is far from providing a realistic picture of 
these sources. A major obstacle in obtaining more complete mod- 
els for the sheared field is non-linearity of the force-free equation. 



Given two force-free fields, Bi and B2, the linear combination 
aBi + bB2 (with a and b two constants) is itself force-free only if 
(V xBi)xB2 + {V xB2)xBi = 0. This implies that a generic 
sheared field can not be expressed as an expansion of sheared mul- 
tipoles, or, conversely, that the superposition of twisted multipoles 
is not a force-free field. An obvious case in which the previous 
condition is satisfied is that of potential fields. Since sheared fields 
depart smoothly from potential multipoles for p ^ po, for small 
enough twists a linear combination of force-free twisted multipoles 
(all with the same twist angle) may provide an approximate force- 
free field. 

Despite many efforts have been devoted to develop techniques 
for solving the force-free equation, V x B = a(x)B, no gen- 
eral, affordable method has been present ed so far. The ca se of a a 
constant has been discussed long ago by IChand raskhar & Kendalll 
(Il957h and more recently bv iMastrano & Melatosi (i2008k in con- 
nection with magnetars. If a is a known function of position, 
ICuperman & Ditkowskil (ll99ll) presented an analytical method for 
solving the force-free equation also in the non-axisymmetric case. 
However, this is of little use for the problem of constructing self- 
consistent force-free magnetospheres since prescribing a is tanta- 
mount to assign the currents which sustain the field, while for the 
case at hand the field and the supporting currents depend on each 
other. A completely g eneral, analytical technique has been pro- 
posed bv lUchidal (1199 7a. b). This is based on a relativistic (tensor) 
description of the electromagnetic field and on the introduction of 
two scalar potentials which are the analogues of the classical Euler 
potentials. It has been show n to be worka ble in the axisymmetric 
case (a non-aligned rotator, IUchidalll998l) and, for the particular 
case of a non-rotating, aligned magnetosphere, it is possible to ver- 
ify that equation ((8]) is recovered. Further work on this is in progress 
and will be reported in a subsequent paper (Pavan et al., in prepa- 
ration). 
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APPENDIX A: 

We start by expressing both the generating function and the eigen- 
value C in terms of a series expansion around the corresponding 
known untwisted quantities (labelled again with the index 0) 

C = Co+(^) Ap + ... 

/(m) =/o(m)+/i(m)Ap + ... 

where Ap = p — po.By substituting the previous expressions into 
equation (IH), we obtain, to first order in Ap, 



(1 - +Po(l +Po)/i(m) + (1 + 2po)/o + (Al) 

It can be easily seen that equation dAlb admits analytical so- 
lutions only for integer values of the exponent 2/po, i.e. for dipoles 
and quadrupoles for which 2/po = 2 and 1, respectively. Since /o 
is itself a solution of the GSS equation, it satisfies the same bound- 
ary conditions we need to impose on /. This implies that the con- 
ditions on /i are, in the case of a dipolar field, /i (1) = f[ (0) = 0, 
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Figure Al. Dipolar angular functions and eigenvalues C(p) obtained 
with fixed flux (upper panels) and fixed Bpoie (lower panels). SoHd lines 
represent the analytical first order approximation, while the dots are the 
numerical solutions of equation ([8). 



supplemented by either /i (1) = or /i (0) = 0. The two solutions 
that obey the previous two sets of conditions are 



and 



dC 
dp 



dC 
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fo- 
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respectively, where /o(m) — 1 — M^- The complete expressions for 
the two generating functions are then 

5/0 + 17 2, 



fw95 /o 1 + 



32 



-fi Ap 



17 
32 



/oAp; 



since Ap < 0, it is fw9^ < /tlk for any value of the parame- 
ter in accordance with numerical results. A comparison of the first 
order approximations with the exact numerical solutions is shown 
in Fig. lAll We find that the agreement between the two is satisfac- 
tory (relative error < 8%) up to A at s = 0.2 in the case of fixed 
intensity and A^ats = 0.4 in the case of constant flux. The ratio 
A = fTLK / fw95 introduced in ^2.2l can be expanded as 



32 ^ 



By applying the same procedure, one can derive the analytic 
first order expansion for quadrupolar fields. For /i(l) = /i(0) = 
/{ (1) = 0, the first two terms in the expansion of / turn out to be 

2\ 



/o(m) 



f [3(1 
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together with 
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